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1  .  INTRODUCTION 

;  Jhe  analytical  study  of  the  nature  of  interlaminar  stresses  near 

a 

stressf-free  edges  has  been  the  subject  of  substantial  research  interest. 
Such  studies  are  significant  because  the  high  interlaminar  stresses  or 
stress  singularities  near  stress-free  edges  may  cause  delamination  failure 
as  shown  by  experimental  invest! gations^Ti'^-  Also  the  accurate  predic¬ 
tion  of  these  stresses  may  be  useful  in  the  design  of  test  specimens  for 
investigation  of  laminate  strength^!.  One  of  the  first  analytical 
studies  on  interlaminar  stresses  near  stress-free  edges  was  published  by 
Pipes  and  Pagano  f5].  In  their  study,  finite  difference  technique  was 
applied  to  a  symmetric  finite-width  composite  laminate  subjected  to  uniform 
uniaxial  strain.  Subsequently  other  techniques  were  used  to  solve  similar 
problems,. Among  these  are  the  boundary  layer  or  the  perturbation  method 
[6-8],  Galerkin  method  X9]  and  the  finite  element  method  tlO-13].  However, 
in  these  studies  the  exact  nature  of  stress  singularities  at  the  free  edges 
were  not  taken  into  account  in  the  formulation.  Therefore,  perhaps  except 
for  Reference  13,  the  accuracy  of  solutions  near  the  free-edge  interface 
appears  marginal  at  best.  In  the  finite  element  analysis  reported  in 
Reference  13,  extremely  fine  meshes  were  used  near  the  free-edge  interface, 
resulting  in  a  finite  element  model  with  a  large  number  of  degrees  of  free¬ 
dom.  Recently  the  exact  nature  of  stress  singularities  near  the  free-edge 
interfaces  in  cross  ply  and  angle  ply  composites  has  been  investigated 
[14-18].  Reference  18  also  includes  interlaminar  stress  distribution 
determined  by  the  boundary  collocation  method. 

In  this  report,  we  present  an  efficient  hybrid  finite  element  method 
for  analysis  of  interlaminar  stress  or  free  edge  effect  in  symmetric 
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composite  laminates  subject  to  uniform  uniaxial  strain.  Both  cross  ply 
and  angle  ply  laminates  are  considered.  The  main  feature  of  the  present 
study  is  the  use  of  a  special  hybrid  element  with  an  embedded  stress  singu¬ 
larity.  The  effectiveness  of  the  hybrid  finite  element  method,  when  applied 
to  problems  with  a  singularity,  such  as  a  crack,  has  been  amply  demonstrated 
by  Pian  and  other  people  [19-21].  Especially  the  bi-material  crack 
analysis  by  Lin  and  Mar  [22]  is  quite  relevant  to  the  present  interlaminar 
stress  analysis. 

In  the  next  section,  a  formulation  for  both  special  and  regular  hybrid 
elements  are  given.  Numerical  examples  are  treated  in  the  third  section. 

A  detailed  discussion  on  the  determination  of  assumed  stress  and  displace¬ 
ment  field  is  included  in  the  Appendix. 
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2.  FORMULATION 


2.1  Problem  Description 

Figure  1  shows  a  long,  symmetric  cross  ply  or  angle  ply  composite 
laminate  loaded  in  the  z  direction.  The  laminate  has  four  plies,  each 
with  thickness  h.  The  width  of  the  laminate  is  2b.  For  the  region 
away  from  the  ends,  the  laminate  can  be  considered  to  be  subject  to  uniform 
uniaxial  strain  [5,23],  and  the  displacement  u,  v  and  w  in 

X,  y  and  z  directions  respectively  can  be  written  as 

u  =  U(x,y) 

V  =  V(x.y)  (1) 

w  =  £gZ  +  W(x,y) 


where  U,  V,  W  are  functions  independent  of  z.  For  cross  ply  laminates 
W(x,y)  *  0.  Also  stress  is  independent  of  z,  and  thus  equilibrium  eqs.  are 


3a 

_ )0(  ^ 

3X 


0 


3a  3a 

0 

3X  3y 


3a  3a 

-^  +  — 0 
3x  3y 


(2) 


Thus  the  problem  reduces  essentially  to  a  two-dimensional  boundary  value 
problem  in  x,y  plane.  For  cross  ply  laminates  a  =  a  =  0. 

y  fc 
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2.2  Hybrid  Finite  Element  Formulation 

In  hybrid  finite  element  approximation,  a  cross  section  of  the  laminate- 
normal  to  the  2  coordinate  is  divided  into  two  regions,  (see  fig.  2). 

The  cross-hatched  region  in  fig.  2  contains  both  the  ply  interface  and  the 
stress-free  edge.  Since  stress  singularity  is  present  at  the  junction  of 
the  ply  interface  and  the  stress-free  edge,  this  region  is  modeled  by  a 
single  special  hybrid  element  with  proper  stress  singularity  terms.  The 
rest  of  the  plane  is  modeled  with  many  ordinary  hybrid  elements. 

For  actual  formulation  of  hybrid  finite  elements,  we  may  start  from 
the  Hellinger-Reissner  principle  or  the  modified  complementary  energy 
principle  [24-25],  Both  are  two-field  variational  principle  with  displace¬ 
ment  and  stress  components  as  independent  variables. 

For  the  Hell inger-Reissner  principle,  the  functional  is  expressed 
as  follows: 

"R  "  1  “ij  •  I  ?  ^ijkt  "ij  “kt  •  j  ’’’i  “i 

V  V  s 

0 

stress  tensor 

i  (u.  .  +  u.  .)  =  strain  tensor  in  terms  of  displacement 

C.  1  >  J  J  >  ' 

vector 

displacement  vector 
compliance  tensor 
applied  traction  vector 
volume  of  the  solid  body 

portion  of  the  surface  where  traction  is  applied 


where 


^ijkx, 

^i  = 
V  = 

s  = 

a 
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The  stationari ty  of  leads  to 


where 


0 

a 


(4) 


aei- 


=  T  (6U, 


1  .J 


(5) 


In  eq.  (4),  6e . .  and  6a..  can  be  interpreted  as  virtual  strain  tensor 
and  virtual  stress  tensor  respectively.  For  the  problem  to  be  considered 
here,  the  traction  vector  =  0.  Now  if  V.|  represents  the  volume  of 
the  special  element  and  ^2  rest,  then  6ttj^  =  0  can  be  expressed  as 


«”R  ’  "r. 


+  6lTr,  =  0 


(6) 


where 


6lTr 


4 


IJ 


6e.  .dV 
IJ 


e .  . ) 
IJ 


6a .  .dV 


(I<  =  1,2) 


(7) 


(A)  Regular  Element 

The  region  away  from  the  free-edge  ply  interface  is  modeled  by  regular 

hybrid  stress  elements  since  no  singularity  is  expected  there.  The  finite 

element  modeling  is  accomplished  by  using  the  expression  for  6t:-  in 

K2 

eq.  (7).  For  the  problem  considered  here,  the  volume  integral  in  eq.  (7) 
reduces  to  an  integral  over  area  A2.  Written  in  matrix  form. 


6Trp  » 


I  6e^  a 


dA  + 


60"^  (e  -  e)  dA 


(8) 
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The  superscript  T  represents  the  transpose. 

Since  is  given,  the  strain  vector  e  and  the  stress 

vector  a  are  expressed  in  component  form  as  follows: 


c 


E 


yy 


(9) 


0 


XX 


yy 


yz 


zx 


xy 


(10) 


For  cross  oly  laminates,  e  =  e  =  a  =  a  =0.  Similar  expressions 

yz  zx  zx 

hold  true  for  6e  and  (5a.  Strain  vector  e  is  related  to  stress  vector 
0  through  the  following  eq. 


E  =  B  a  +  E„ 
-  -  ~o 


(11) 


where  B  is  now  a  modified  compliance  matrix,  and  e^  plays  the  role  of 
an  initial  strain  vector.  See  Appendix  A.l  for  derivation  of  B  and  e^. 
Substituting  eq.  (11)  into  eq.  (8) 


a  dA  + 


B  a  - 


(12) 
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For*  finite  element  approximation,  the  displacement  vector  u  is  assumed  in 
each  element  as 


y  =  N  (13) 

where 

N  =  shape  function  matrix 
=  element  nodal  displacement  vector 

Then  symbolically 
?  =  ?  Se 

and  also 

«J  =  f5ge 

The  assumed  stress  field  satisfies  equilibrium  in  each  element,  and  can 
be  expressed  as 

0  =  P  e  (16) 

where 

P  =  stress  shape  function  matrix 
B  =  the  vector  of  unknown  stress  parameters 

Also, 


=  P  6  B  (17) 

Substituting  eqs.  (13)  -  (17)  into  eq.  (12) 

6it„  »  Z  [6q„^  6  +  66*'^  (G  q„  -  H  B  -  G„)]  (18) 
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where 


G  =  j  p''’  f  dA 

H  =  I  p"'’  B  P  dA  (19) 

=  1  P"’’  dA 

-0  J  '  -0 

The  I  notation  indicates  summation  over  all  regular  elements.  From 
eq.  (18), 

e  Qe  -  H  B  -  Gq  =  0  (20) 

for  arbitrary  66.  From  eq.  (20) 

5  =  «■'  (5  9e  ■  9o> 

Substituting  eq.  (21)  into  eq.  (18) 

^  9e  '  §0^ 

(22) 

=  ^  (kg  9e  ■ 

where 

k  =  g'*’  G  (23) 

~e  ~  ' 

is  the  element  stiffness  matrix  and 

Q  =  G„  (24) 

3e  -  -  '0 
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is  the  element  nodal  load  vector  due  to  prescribed  strain 

The  regular  elements  in  the  present  study  are  four  node  element  with 
isoparametric  representation  for  the  assumed  displacement  field.  The 
assumed  stress  field  is  linear  and  satisfies  equilibrium  equation  within  - 
each  element.  Thus  for  the  cross  ply  cases,  a  regular  element  has  eight 
nodal  degrees  of  freedom  and  the  following  stress  field  with  seven  stress 
parameters : 


“xx '  8]  *  Sa*" 


“yy  °  *  “s’"  *6^ 

“xy  '  ®7  -  ®6’'  - 


For  angle  ply  case,  an  element  has  twelve  nodal  degrees  of  freedom.  For 
the  assumed  stress  field,  the  following  and  components  are 

added  to  those  in  eq.  (25) 


xz 


69X 


"^10^ 


yz 


=  6 


11 


612X 


'  BgY 


(26) 


(B)  Special  Element 

The  special  element  incorporates  stress  singularity.  In  addition,  the 
stress-free  condition  along  the  edge  and  the  bonding  condition  along  the  ply 
interface  are  exactly  satisfied,  (See  fig.  4.). 

Using  the  divergence  theorem,  Sup  in  eq.  (7)  can  be  transformed  to 


9 


r  _  ^ 

T.  6U.  + 

h 


6T.  (u.  -  u.)  dS 


‘1 


(27) 


for  stress  which  is  in  equilibrium  and  also  compatible.  In  eq.  (27) 

is  the  displacement  integrated  from  stress  and  is  independent  of  . 
Written  in  matrix  form. 


611 


T  dS 


(u  -  u)  dS 


Since 


and 


u 

<•  \ 

U(x,y) 

U  =  ' 

V 

'  =  i 

V(x,y) 

.  w  . 

CqZ  +  W(x,y) 

u 

✓  N 

U(x,y) 

u  =  ' 

7 

=  < 

V(x,y)  7 

.  w  .. 

t  z  +  W(x,y) 

^  y 

(STTn  can  be  rewritten  as 
*^1 


6irr 


I  ds  +  j  6T^  (y  -  y)  ds 

h  h 


(28) 


(29a) 


(29b) 


(30) 
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where 


U  = 


(31) 


and 


f  T  ^ 

X 

°XX 

°xy 

'^X2 

T  =  ^ 

T 

y 

>  = 

^xy 

°yy 

'^yz 

< 

.  ‘^X2 

°yz 

°22. 

<  ^  > 


(32) 


In  eq.  (32),  z,  m  and  n  are  the  components  of  a  unit  vector  normal  to 
the  surface.  Similar  expressions  hold  for  U,  sU"  and  6T.  In  eq.  (28), 
integration  is  defined  over  the  surface.  However,  for  the  present  problem, 
it  reduces  to  a  line  integral  along  the  element  boundary. 

For  finite  element  approximation,  it  is  convenient  to  separate  g  and 
U  such  that 


q  =  *g  + 

U  =  *U  +  °U 


(33) 


Here  term  is  constant  stress  predetermined  to  take  care  of  the 
term  in  eq.  (11).  The  *q  terms  represent  the  assumed  stress  and  are 
free  of  terms.  The  °y  and  *y  vectors  correspond  to  °q  and  *g 
respectively.  The  pair  °U  and  °o  satisfies  equilibrium,  compatibility 
and  stress-free  conditions  along  the  edge  as  well  as  the  bonding  condition 
along  the  ply  interface.  See  Appendix  A. 2  for  details.  The  pair  *y  and 
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*a  also  satisfies  all  these  conditions. 

Symbolically,  the  assumed  stress  *q  can  be  expressed  as 

*q  =  P  I  (34) 

where  6  is  the  vector  of  unknown  parameters.  The  *q  vector  includes 
singular  as  well  as  regular  stress  terms.  From  eqs.  (32)  and  (33),  we  may 
write 

T  =  *T  +  °T  (35) 

And  then,  from  eq.  (32)  and  (34),  *T  may  be  written  symbolically  as 


*T  =  R  e  (36) 

The  *y  vector  is  also  expressed  symbolically  as 

*y  =  L  8  (37) 


The  displacement  vector  U  is  assumed  in  terms  of  element  nodal  displace¬ 
ment  vector  q^  such  that 

U  =  N  (38) 


The  N  matrix  now  represents  the  shape  function  matrix  along  the  element 
boundary. 

With  eqs.  (36)  to  (38),  Sttd  in  eq.  (28)  can  be  written  as 

*^1 


"  ^5e^  5o  *  9e  ■  -5  ■  9o^ 


where 
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N  dS 


L  dS 


R^  °X  dS 


R^  °y  dS 


(40) 


The  H  matrix  is  symmetric  although  the  integrand  R^  L  is  not.  It 
should  be  noted  that  the  path  of  line  integral  in  eq,  (40)  does  not 
include  the  stress-free  edge  and  the  ply  interface.  Thus  the  integration 
path  does  not  cross  the  singular  point.  In  fact,  this  is  the  essence  of 
hybrid  formulation  for  analysis  of  cracked  solids  [19].  From  eq.  (39) 

G  H  6  -  =  0 

-  ^e  -  ^  -0 

or 

§  = 

Substituting  eq.  (41)  into  eq.  (39), 

5^  «■’  (5  9e  -  ?o 

(42) 

°  (5e  3e  -  ?e> 
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where 


kg  "  §  ^^3) 

is  the  stiffness  matrix  of  the  special  element  and 

?e  '  5^^  t!"'  5o  ■  9o 

is  the  element  load  vector. 

The  special  element  used  in  the  present  study  has  nine  nodes  as 
shown  in  Fig.  3.  The  displacement  between  two  nodes  is  linear  and  thus 
compatible  with  the  adjacent  regular  elements.  The  number  of  unknown 
parameters  in  the  0  vector  is  three  for  cross  ply  case  and  four  for 

angle  ply  case.  See  Appendix  A. 3  for  details  on  the  assumed  stress  field. 

(C)  Summing  or  Assemblying 

The  finite  element  equation  for  the  whole  problem  is  derived  by 
summing  airn  and  6irD  such  that 

Ki  K2 

Now  the  summation  notation  stands  for  summing  or  assemblying  over  all 
elements  and 

g  =  global  nodal  displacement  vector 
k  =  global  stiffness  matrix 

Q  *  global  load  vector 
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For  arbitrary  6q,  we  obtain 


k  q  =  Q 

which  can  be  solved  for  q. 
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3.  NUMERICAL  EXAMPLES 

The  effectiveness  of  the  present  method  has  been  tested  by  solving 
three  examples  of  four-ply  laminates.  They  are  [90°/0°]g  and  [0®/90'']g 
cross  ply  laminates  and  a  [±45°]^  angle  ply  laminate  as  shown  in  fig.  5. 

The  geometrical  dimensions  and  material  properties  used  in  the  present 
study  are  given  as  follows: 

(a)  geometry 

half  width  b  =  24" 
ply  thickness  h  =  3" 

(b)  material  property 
E^^  =  20.0  X  10®  psi 

E22  “  2.1  X  10®  psi 

''12  '  ''23  "  '’31  ' 

^1 2  *  ^23  ~  ^31  ~  ^  ^ 

These  properties  are  the  same  as  those  used  by  other  people  [5,9,12].  Due 
to  symmetry,  only  a  quarter  of  the  section  (the  upper  left  part)  was  modeled. 
Two  different  meshes,  coarse  and  fine,  were  used  to  check  convergence. 

Figure  6  shows  these  meshes  in  scale.  The  number  of  nodes  for  the  coarse 
and  fine  meshes  are  111  and  159  respectively. 

(A)  Cross  Ply  Case 

Figures  7  and  8  shows  stress  distribution  along  the  90“-0®  ply 
Interface  at  y  *  h  for  both  [90®/0‘*]^  and  [0“/90‘’]^  laminates.  These 
results  are  for  the  fine  mesh.  Although  they  are  not  shown,  solutions 
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obtained  by  the  coarse  mesh  are  very  close  to  those  by  the  fine  mesh, 
indicating  convergence.  Note  that  for  convenience  a  new  coordinate  x  is 
introduced  as  shown  in  fig.  2.  The  x  coordinate  is  introduced  such  that 
X  -  0  at  the  free  edge  and  x  =  b  at  the  center. 

Another  measure  of  convergence  is  to  check  the  value  of  free  parameter 
S-j  corresponding  to  the  singular  stress  term.  For  the  [90°/0°]g  laminate, 

Q 

the  computed  values  of  6^  are  0.1044  x  10  for  the  coarse  mesh  and 
0.1036  X  10®  for  the  fine  mesh.  For  the  [0°/90°]^  laminate,  these 
values  are  -  0.4384  x  10®  for  the  coarse  mesh  and  -  0.4338  x  10® 
for  the  fine  mesh. 

In  fig.  7,  normal  stress  shows  high  gradient  near  the  free  edge 

for  both  [0®/90]  and  [90°/0°]  laminates.  The  [90°/0°]  case  is 
s  s  s 

particularly  interesting.  Results  reported  by  Wang  and  Crossman  [9]  and 
Spilker  [12]  indicates  very  small  at  the  free  edge.  On  the  other 

hand,  the  present  result  shows  ever-increasing  hi  ah  positive  normal  stress 
confined  to  extremely  narrow  region  at  the  edge.  This  result  is  in  agree¬ 
ment  with  that  by  Raju  [13],  and  clearly  indicates  existence  of  stress 
singularity.  It  appears  that,  although  Wang  and  Crossman,  and  Spilker 
used  very  fine  meshes  near  the  edge,  the  size  of  elements  was  not  small 
enough  to  capture  the  detailed  picture.  Meanwhile  the  size  of  elements 
used  in  reference  13  at  the  edge  was  extremely  small,  resulting  in  a  model 

with  an  excessively  large  number  of  unknowns.  Figure  8  shows  a  distri- 

xy 

bution.  It  attains  maximum  value  very  close  to  the  free  edge.  Of  course 
it  drops  to  zero  at  the  free  edge  itself. 

(B)  [±45°],  Angle  Ply  Case 

8 

In  this  example,  computed  values  of  6^  are  0.1131  x  10  for 

Q 

the  coarse  mesh  and  0.1123  x  10  for 


17 


the  fine  mesh,  indicating  convergence.  Figures  9  and  10  show  the  computed 
stress  distribution  along  the  ply  interface  at  y  =  h.  Figure  9  indicates 
ever-decreasing  (negative)  near  the  free  edge,  indicating  singular 

stress.  However,  approaches  zero  as  x  increases.  In  addition,  shear 

stress  shows  an  ever-increasing  trend  toward  the  free  edge.  Other  stress 

components,  o  ,  a  ,  a  are,  of  course,  zero  at  the  free  edge.  Shear 

XX  X2 

stress  almost  constant  away  from  the  free  edge,  in  accordance  with 

the  classical  lamination  theory.  It  attains  maximum  near  the  free  edge 

before  it  drops  to  zero.  The  normal  stress  a  also  reaches  maximum 

near  the  free  edge  and  then  quickly  reduces  to  zero  at  the  free  edge  itself. 

The  result  for  shear  stress  a  is  not  presented  here  because  its  magni- 

xy 

tude  is  very  small  compared  with  other  components. 

Figure  11  shows  a  along  the  free  edge.  Here  we  observe  a  sharp 

•J  V 

change  in  magnitude  near  the  ply  interface.  This  behavior  is  consistent 
with  the  existence  of  stress  singularity.  The  results  presented  here  agree 
with  those  in  reference  13.  However,  in  reference  13,  an  extremely  fine 
mesh  had  to  be  used,  and  the  exact  nature  of  stress  singularity  could  not 
be  determined  within  reasonable  accuracy. 


4.  CONCLUSION 


Numerical  results  indicate  that  the  hybrid  finite  element  formulation 
involving  a  special  element  with  embedded  stress  singularity  is  a  very 
efficient  means  for  accurate  determination  of  interlaminar  stress  distri¬ 
bution.  For  both  cross  ply  and  angle  ply  symmetric  laminates  considered 
here,  the  present  method  provides  converged  stress  values  near  the  junction 
of  the  stress-free  edge  and  the  ply  interface.  These  stress  values  are 
much  more  accurate  than  those  obtained  by  others  using  conventional  finite 
element  models  that  do  not  include  proper  singularity.  With  the  present 
formulation,  it  is  possible  to  use  a  much  more  coarse  finite  element  mesh, 
resulting  in  a  substantial  improvement  in  computing  effort. 
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Finite  elemen 


APPENDIX 


A.l  ANISOTROPIC  ELASTICITY  WITH  PRESCRIBED  UNIFORM  UNIAXIAL  STRAIN 

For  three-dimensional  solids,  the  three  equations  of  equilibrium  are 

expressed  in  terms  of  stress  components  a  ,  a  ...  a  ; 

XX  yy  xy 


9cr.,..  3a...,  3a 


XX 


+  0 


3X  3y  3Z 


3a  3a  3o 
3X  3y  32 


(A. 1.1) 


!!xi  +  +  , 0 

3X  3y  32 


The  six  strain  components  e  ,  e  ,  .  .  .  e  are  related  to  the 

xa  yy  xy 


six  stress  components  as  follows: 


^xx  *  ^11  '^xx  ^12  V  ^13  ^22  ^14  V  ^15  '^zx  ^16  ""xy 
Vy  "  ^21  ^xx  *  ^22  ^yy  *  ^23  ‘^zz  ^24  ^25  '^zx  ^26  °xy 


^zz  "  ^31  °xx  ^  ^32  °yy  ^33  "^zz  ^34  ^yz  ^35  °zx  ^36  °xy 
^yz  =  ^41  °xx  ^  ^42  V  ^  S43.'’zz  ^  ^44  V  ^  ^45  ^zx  ^  ^46  ^xy 


(A. 1.2) 


e  -  Sri  a  +  Sco  a  +  Sr-,  a  +  Sc>,  a  _  +  Srr  a  +  Srr  a 

2X  51  XX  52  yy  53  zz  54  yz  55  zx  56  xy 


e  ~  Sri  a  +  Sro  a  +  Sr-,  a  +  Sr*  a  +  Srr  a  +  Srr  a 

xy  61  XX  62  yy  63  zz  64  yz  65  zx  66  xy 


34 


where  S^2  ®tc.  are  compliance  coefficients.  The  strain-displacement 

relation  is  expressed  as 


3u 

XX 

3X 

3V 

yy 

ay 

3w 

22 

32 

3V 

3W 

yz 

3Z 

ay 

Iw 

4. 

3U 

2X 

3X 

32 

3U 

3V 

xy 

ay 

3X 

(A. 1.3) 


where  u,  v  and  w  are  displacement  components  in  x,  y  and  z  directions 
respectively. 

If  the  body  is  subjected  to  uniform  uniaxial  strain  =  Eq,  stress 
components  do  not  vary  along  the  z  direction  [23].  Then  the  equilibrium 
equation  reduces  to 


3o 

3a 

XX 

+ 

3X 

ay 

3a 

3a 

_XZ 

+ 

-TL 

3X 

ay 

3o 

3o 

X2 

+ 

3X 

ay 

Now  from  the  third  equation  of  (A. 1.2) 


'"zz  “  ^  ■  ^31  °xx  “  ^32  °yy  '  ^34  V  '  ^35  '’zx  ‘  ^36  °xy' 

(A. 1.5) 
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Substituting  (A. 1.5)  into  the  remaining  five  equations  in  eq.  (A. 1.2),  we 
obtain 


'xx  =  ®ll“xx  ®12“yy  *  *  hs^zx  “le’xy  *  ^0 


=yy  °  ^ai^xx  *  ®22‘’yy  ®24''yz  *  ®25“2x  *  ®26°xy  *  S33  ‘o 


:  =  B/nO  +  B.oO  B/i/,a  _  +  B*rO_„  +  B/,cO  ■*■  c —  (A.1  .6) 

yz  41  XX  42  yy  44  yz  45  zx  46  xy  0 


^5  3 

^zx  "  ®51°xx  ®52°yy  *^54^  ®55°zx  ^56°xy  ^  ^0 


^xy  "  ®61^xx  ®62''yy  ^64°yz  ^65®xz  ®66°xy  S33  ^0 


where 


S.,  S., 

B.  .  =  S.  - 


ij  ij  S 


(i,j  =  1,2, 4, 5, 6) 


33 


(A. 1.7) 


For  symmetric  angle  ply  case 


Bi4  =  B24  =  B34  =  B54  =  0 


®16  ^  ®26  ^  ®36  ""  ^'6  "  ° 


In  matrix  form,  eq.  (A. 1.6)  can  be  written  as 


(A. 1.8) 


e  =  B  o  +  e, 


(A. 1.9) 
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wi  th 


e  - 


^  t  '' 
XX 


yy 

^  'yz  ^ 


zx 


xy 


B  = 


Ml 


’61 


=  < 


XX 


yy 


yz 


zx 


'^xy 


^  ^13  ^ 


So  ^  1 


’23 


’43 


’53 


'63 


(A.l .10) 


’16 


’66 


=5x5  matrix 


(A. 1.11) 


(A. 1.12) 


(A. 1.13) 


’33 
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For  the  present  problem,  the  displacement  components  u,  v,  w  can 
be  expressed  as 

u  =  U(x,y) 

V  =  V(x,y)  (A.l .14) 

w  =  CqZ  +  W{x,y) 


Then  the  strain-displacement  relation  is  written  as 


=  m 

^xx  ”  ax 


-  li 
^yy  '  3y 

_  ^ 
^yz  ■  ay 


- 

^zx  ”  ax 


=  iU  +  iV 

^xy  ”  3y  ax 


(A.l  .15) 


For  cross  ply  laminate,  further  simplification  is  possible  since 


a  =  a  =0 
yz  zx 


G  =  e  =0 

yz  zx 


(A.l  .16) 
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The  equilibrium  equation  now  reduces  to 


do 


XX 

3X 


0 


3a  3a 

+  -M. 


3X 


sy 


0 


Also 

^13 

^xx  "  °xx  ®12  ""yy  *  ^0 

^23 

^yy  "  ®21  ""xx  ®22  °yy  ^0 


Sy  '  ®66  '"xy 
and 


^  3U 

^XX  ~  3X 

=  — 

^yy  ~  3y 

=  ly.  +  iv 
^xy  "  3y  3x 


(A. 1.17) 


(A. 1.18) 


(A. 1.19) 
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A. 2  CONSTANT  STRESS  FIELD** 

The  role  of  predetermined  constant  stress  field  is  to  separate  terms 
involving  prescribed  strain  from  the  assumed  stress  field  and  the 
corresponding  assumed  displacement  field.  For  the  present  special  element, 
the  constant  stress  term  satisfies  the  stress-free  condition  along  the  free 
edge.  In  addition,  the  constant  stress  field  and  the  displacement  field 
integrated  from  the  constant  stress  field  satisfy  the  bonding  condition 
along  the  ply  interface. 

A. 2.1  Cross-Ply  Case 

For  the  coordinate  system  shown  in  figure  4,  the  stress-free  condi¬ 
tions  are; 


at  ®  i 

and 


at  ®  "  ■  f 

The  bonding  condition  along  6=0®  are 

u  _  i 

a  ~  o 
xy  xy 

u  _  l 

°yy  ~  ^yy 


(A. 2. 2) 


(A. 2. 3) 


**  The  senior  author  has  recently  noticed  a  similar  development  in  ref.  17. 
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The  superscripts  u  and  £  stand  for  the  upper  layer  and  the  lower  layer 

respectively.  The  constant  stress  components  ,  °a  ,  °a  that 

XX  yy  xy 

satisfy  the  above  stress-free  and  the  bonding  conditions  are 


0  u  0  £ 

a  =  a 
XX  yy 


=  0 


0  U  0  £ 

a  =  a 
xy  xy 


=  0 


(A. 2. 4) 


0  U  0  s, 
a  =  a 

yy  yy 


=  c 


where  c  is  a  constant  to  be  determined  as  follows.  Substituting  the 
constant  stress  terms  in  eq.  (A. 2. 4)  inco  eq.  (A. 1.18)  and  integrating, 
we  obtain  the  corresponding  displacement  components  and  for 
the  upper  layer  as  follows 


V  =  (B,2“  c  X 

^33 


°v'‘  "  (822“  c  t  ^  t„)  y 


U  0 ' 


^33 


(A. 2. 5) 


excluding  the  rigid  body  modes.  Similar  expression  holds  for 

0  0  £ 

displacement  components  U  and  V  .  Then  from  the  displacement 
continuity  at  the  90°  -  0°  ply  interface. 

S  S  ^ 

(612“  c  +  -^  Cq)  X  =  c  +  Eq)  X  (A. 2. 6) 

^33  ^33 
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Solving  for  c,  we  obtain 


c 


'13 


'33 


(A. 2. 7) 


A. 2. 2  Angle-Ply  Case 

Now  stress-free  conditions  are: 


u  u 
0  =0 
XX  xy 


(A. 2. 8) 


at  ^  ^  I 

and 


a 


(A. 2. 9) 


at 


e  = 


TT 

2 


The  bonding  conditions  along  the  ply  interface  (e  =  0°)  are; 


c 


U  I 

a  =0 

yy  yy 


and 


(A. 2. 10) 


(A. 2. 11) 
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The  constant  stress  components  that  satisfy  the  above  stress-free  and  the 
bonding  conditions  are 


(A. 2. 12) 


Substituting  this  stress  state  into  eq.  (A. 1.15)  and  integrating,  we  obtain 
the  corresponding  displacement  fields.  Then  from  the  matching  condition 
eq.  (A. 2. 11)  of  displacements  along  the  ply  interface 


(A. 2. 13) 


and 


0  U  0  I! 
a  =0 

yy  yy 


(B52  -  Bp,  )  S,,u 


c  S  ^ 


(A. 2. 14) 


52 


'33^ 


’53 


or  since 


(A. 2. 15) 
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for  synnnetric  angle  ply 


0  u  _  0  j  _  ^0  ^53 

0....  -  0  -  ~  - 


yy  -yy  r  S  “ 
“52  “33 


In  addition. 


0||U  _/dUO^  u  13  ^  )  y 

u  -  (b^2  V  ^ 

“33 


o.,u  _  /n  u  o„  u  23  ^  \  „ 

^  -  ^^22  V  y 

“33 


°w“ 


=  0 


Similar  expressions  hold  for  stress  '  '  '  ^'’xy^ 

°U^,  °V^  and  V  for  the  lower  layer. 


(A. 2. 16) 


(A. 2. 17) 


displacement 


A. 3  ASSUMED  STRESS  FIELD  (SINGULAR  AND  NONSINGULAR) 

The  assumed  stress  and  the  corresponding  displacement  fields  satisfy 
all  governing  equations  of  elasticity,  the  stress-free  condition  over  the 
free  edge  and  the  bonding  condition  along  the  ply  interface.  In  order  to 
ensure  convergence,  the  assumed  stress  must  include  singular  terms  [26]. 
The  determination  of  singular  stress  fields  proceeds  as  follows. 


A. 3.1  Cross-Ply  Case 

By  introducing  a  stress  function  F(x,y)  such  that 


XX 


yy 


xy 


lif. 

3y^ 

3^F 


(A. 3.1) 


3X 


jIl 

3x3y 


the  equilibrium  equation  (A. 1.17)  is  always  satisfied,  and  the  compatibility 
among  strain  components  in  eq.  (A. 1.18)  leads  to  the  following  fourth-order 
equation  for  F. 


822  '^5  "  <2  8,2  *  Bjj)  =  0 


(A. 3. 2) 


The  above  equation  has  the  general  solution  of  the  following  form: 


F(x,y)  =  F|^  (x  +  pj^y) 


(k  =  1.2. 3,4) 


(A. 3. 3) 


where  P|^  is  a  root  of  the  following  fourth  order  algebraic  equation. 
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(A. 3. 4) 


B^1  +  (2  B^2  ^56 )  ^  ^22  ~  ® 

The  solutions  to  the  above  equations  constitute  two  conjugate  pairs.  Thus 
if  and  y2  are  two  distinct  roots,  then 


"3  =  "l 


(A. 3. 5) 


where  is  the  conjugate  of  u^j  etc.  In  order  to  determine  the  singu¬ 
lar  stress  field  and  also  the  non-singular  stress  fields,  we  express 
as 


F^(x  t  „|^y) 


\  (a  +  2)(a  +  l) 


(A. 3. 6) 


where  is  a  coefficient  and  o  is  a  quantity  to  be  determined  by  an 
eigenvalue  analysis.  The  local  coordinates  x  and  y  are  related  to 
the  polar  coordinates  r  and  e  as  follows  (fig.  4). 

X  =  r  cose 

(A. 3. 7) 

y  =  r  cose 

Then 


X  +  y,^y  =  r  C 


(A. 3. 8) 


where 


C|^  =  cose  +  sine 


(A. 3. 9) 
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From  eq.  (A.3.1),  the  non-constant  stress  components  are  expressed  as 


4  2  a 

u  _  a  _  .  u  ,  Uv  u. 

XX  "  ^  \  ^  ^ 


\  ^k  (C,,  ) 


(A. 3. 10) 


for  the  upper  layer.  The  corresponding  displacement  components  are  deter¬ 
mined  by  integrating  eq.  (A. 1.18) 


(A.3.11 ) 


where 


Pk  "  ^11  ^k^  ^  ®12 
'^k  "  ^2  ®22^^k 


(A. 3. 12) 


Similar  expression  holds  for  the  lower  layer.  Note  that  term  does 
not  appear  in  the  above  equation  since  it  was  taken  into  account  by  the 
constant  stress  field  and  the  corresponding  displacement  in  the  previous 
section. 
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Applying  the  bonding  condition  along  0=0, 


4 

r 

k=l 


a 


4 

I 

k=l 


(c/) 
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Z 

k=l 


a 


4 

I 

k=l 


(A. 3. 13) 
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k=l 


(C  ") 


a+1 
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Z 

k=l 


(C/) 


a+1 


The  above  equation  can  be  rewritten  in  matrix  form  as 


where 


/ 


A 


u 


> 


(A. 3. 14) 


(A. 3. 15) 
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(A. 3. 16) 


and  T  is  a  4x4  matrix, 
e  =  |-  (cose  =  0,  sine  =  1 ) , 


Since  ^  =  *o  =  0  at 
XX  xy 


4 

z 

k=l 


2+a 


0 


4 

Z 

k=l 


(A. 3. 17) 


And  also  =  *o  ^  =  0  at  9  =  -y  (cose  =  0,  sine  =  -1) 

X  X  xy  w 


4 

z 

k=l 


2+a 

) 


=  0 


4 

Z 

k=l 


(A. 3. 18) 


Equations  (A. 3. 17)  and  (A. 3. 18),  in  conjunction  with  eq.  (A, 3. 14), leads  to 
a  system  of  homogeneous  equations  for  A*^  as: 
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C(a)  =  0 


(A. 3. 19) 


where  C(a)  is  4x4  matrix  with  the  unknown  a.  For  nontrivial  A^, 
the  determinant  of  C(a)  must  be  equal  to  zero,  which  leads  to  an  eigen¬ 
value  problem.  The  eigenvalue  a  can  be  either  real  or  complex.  For 
Singular  stress, 

-1  <  Re(a)  <  0  (A. 3. 20) 

where  Re(a)  is  the  real  part  of  a.  For  the  material  properties  used  in 
the  present  problem,  the  first  eigenvalue  corresponding  to  stress  singularity 
is  a  =  -  C. 333888. 

Us’ng  eigenvalues  and  eigenvectors  determined  from  eigenvalue  analysis, 
we  may  construct  the  assumed  stress  and  displacement  field  as  follows: 


"  a .  ** 

E  r  ^  Re  [  I  (A  ^) 
i=l  k=l 


if 


0 


a .  4 

Z  r  1  Re  [  I  (a/")  (C, '')  ]  - . 

=  1  k  =  l  "  a=a.  ^ 


(A. 3.21  ) 


*0 
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xy 


i=l 


4 

Re  [  l 
k=l 


(A^")  _ 

K  a“a,j 
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and 


^  ^  K'>  h" 

1  =  1  1  k=l  a=ci . 


)  ]  6. 


ct  *+1 

n 


(A. 3. 22) 


1-1  1  k=l  a=a . 


where  ,  ...  6^  are  unknown  free  parameters.  In  the  above  expression, 
the  following  relations  are  utilized. 


^k+2  ^  ^k 


(k  =  1,2) 


(A. 3. 23) 


and  thus 


Pk+2  '  \ 


^k+2  '  ^k  ^ 


(A. 3. 24) 


^k+2  ^  ^k 


In  the  present  study,  a  three  term  (n  =  3)  approximation  was  used. 


A. 3. 2  Angle  Ply  Case 

The  stress  components  are  now  expressed  in  terms  of  stress  functions 
F(x,y)  and  \j;(x,y)  such  that 


(cont ’d) 
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xy  9x3y 

3iJj 

^xz  “  3y 

_  _  ^ 
^yz  “  '  3x 


(A. 3. 25) 


The  compatibility  among  strain  components  in  eq.  (A. 1.5)  leads  to  the 
following  two  homogenous  equations  for  F  and  ij;: 


L.F  +  L,!j/  =  0 


(A. 3. 26a) 


L3F  +  124*  =  0 


(A. 3. 26b) 


where  the  differential  operators  L2»  and  are  given  as 


U  =  B, 


2B, 


+  B. 


■2  “44  ““45  3X3y  “55  ^^2 


3^  3^  33  3 

"  ®24  3  ^  ^^25  ^  ^46^  2  ^^14  ^  ^56^  ?  ®1 5  3 

^  3x  3x  3y  '  3x3y‘^  ‘ 


*-4  ^^22  7^"  7X~  (2B.J2  +  B66)  2  2 

3X  3X  3y  3x  3y 


3y 

(A. 3. 27) 


Combining  eqs.  (A. 3. 26a)  and  (A. 3. 26b)  leads  to 
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(L44  -  Lj^)  F  =  0 

The  above  equation  has  the  general  solution  of  the  form 
F(x,y)  =  (x  +  u^y)  (k  =  1.2, ...6) 

where  is  a  root  of  the  following  sixth  order  equation; 
t^{u)  ^2(11)  -  =  0 

with 

2.2(u)  =  u  -  u  + 

3  2 

t3(y)  =  y  -  (B^4  +  Bgg)  y  +  (B25  +  B4g)y  -  B24 
^4(v)  *  B^^  y^  -  2B^g  y^  +  (28^2"^  '  ^^26  ^  Hi 


The  solutions  to  the  eq.  (A. 3. 30)  constitute  three  conjugate  pairs 
fore  if  y-j ,  y2  and  y^  a<"e  three  distinct  roots,  then 


^'k+3  " 


(k  =  1,2,3) 


Also  from  eq.  (A. 3. 26b)  the  solution  for  can  be  expressed  as 

d  F|^(x  +  y|^y) 

'^k^^  ^k^^  "  “  ^k  d{x  +  yj^y) 


where 

i  =  _  ^3^^k^ 

k  Ji2^^|(^ 


(A. 3. 28) 


(A. 3.29) 


(A. 3.30) 


(A. 3. 31) 


The re - 


(A. 3. 32) 


(A. 3. 33) 


(A. 3. 34) 
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Again,  using  polar  coordinates  r  and  e  (fig,  4) 
X  +  y  =  r  (cos  e  +  sine)  =  r 


where 


C|^  =  cos  e  +  sin  e 


Introducing  and  into  eq.  (A. 3. 25) 


6  2  ot 

”xx  -  ''k  (“k  >  <'^k  > 


*o  “  =  r’  I  A|,“  (0.“)° 

yy  ,^=1  k  ^  k  ' 


*a  -  r“  I  A  u  “  (C  ”)“ 


6  a 

*„  U  _  a  -  .  U  .  U  U  /p  Ux 

°xz  "  ^  \  \  ^^k  ^ 


*“yz“  =  j,  ('k“)“ 


The  corresponding  displacements  are  determined  by  integrating  eq. 
without  Eq  terms.  They  are  written  as 


(A. 3. 35) 


(A. 3. 36) 


(A. 3. 37) 


(A. 1.6) 


(cont’d) 
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(A. 3. 38) 


where 


_  r“*'  ?  ,  u  u  ,, 

''  -  ^TT  \  "k  <'^k  > 


a+l  6  ct+l 

*w“  =  r_  z  (C^) 


u  2 

„  _dU/  Ux.n  u  D  u  u.  ,  u,o  u  u  B 
Pk  "  ^^k  ^  ^12  ‘^16  ^k  ^k  ^^15  '"k  ‘^14  ^ 


Boo 

R  ^  II  ^ 

n2  ^k  u 


B, 
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B,  "  +  (B„,‘^  -  -^) 

26  k  '  25  li.  ' 


B  LI  g  u 

u.n  u  u  24  „  u  u  /g  LI  44  X 

•^k  '  n4  '"k  u  ®46  ^k  ^^45  ■  u  ^ 

^k  ^k 


Similar  expressions  hold  for  the  lower  layer  with  coefficients 
From  the  bonding  conditions  along  e  =  0°, 
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(A. 3. 39) 


£ 

k  * 


(cont'd) 
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Z  A. p  “ 
k=l 


”  -  £  £ 
^  \  Pk 
k=l 
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^-u  u_5.£  £ 

k=l  k=l 


The  above  six  equations  provide  a  relationship  between  and 


From  the  stress-free  conditions  at  e  =  t 


6  2+a 

k=l 


I  A  “  (u  “)  =  0 

k=l  ^ 


1+a 


6  1+a 

V  ^  =  ° 


and  also  from  the  stress-free  conditions  at  0  =  ■  f 


6  .  .  2+a 

I  a/  (-y/)  =  0 

k=l 


t  A  ‘  (-„  *)  ■  0 

k»l 


1+a 


(A. 3. 40) 


(A. 3. 41) 


(cont'd) 
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(A. 3. 42) 


6 

I 

k=l 


Combining  eqs.  (A. 3. 40),  (A. 3. 41)  and  (A, 3. 42),  we  obtain  a  system  of 
homogeneous  equations 

C(a)  A^^  =  0  (A. 3.43) 


where 


A 


u 


(A. 3. 44) 


and  C(a)  is  a  6x6  matrix  with  eigenvalue  a.  For  the  material  proper¬ 
ties  used  in  -his  study  and  ±45®  layer,  the  first  eigenvalue  corresponding 
to  stress  singularity  is  a  -  -  0.0255756. 

Using  eigenvalues  and  the  corresponding  eigenvectors,  we  may  construct 
the  assumed  stress  and  displacement  fields.  For  example,  for  assumed 
stress 

n  a .  6  *i 

=  E  Re  [r  ^  Z  (A^^)  (C,^‘^)  ]  6.  (A. 3. 45) 

i  =1  k=l  a=a . 


wi  th 


2 


The  other  assumed  stress  components  have  a  similar  expression  with  proper 
definition  of  as  follows: 
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Stress  Components 

*  u 
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For  assumed  displacement  *U“, 


II  ^  1  ® 

*U^  =  E  Re  E 

i=l  “r'  k=l 


(V) 


a=o 


U  / U\ 

Pk  e^k  > 


]  B, 


(A. 3. 46) 


7 


The  expression  for  *V^  is  obtained  by  replacing  Pj^*^  in  eq.  (A. 3. 46) 
with  Likewise,  the  expression  for  *w'^  is  obtained  by  replacing 

P|^^  with  rj^^.  Similar  expressions  hold  for  ...  *U^,  *V^  and 

fl 

*W  for  the  lower  layer.  The  coefficients  61,62.  .  .  ST'S  unknown 
free  parameters.  In  the  present  study  a  four  term  approximation  (n=4) 
was  used. 
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